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Abstract 

Two mathematical models of macroevolution are studied. These models have population dynam- 
ics at the species level, and mutations and extinction of species are also included. The population 
dynamics are updated by difference equations with stochastic noise terms that characterize pop- 
ulation fluctuations. The effects of the stochastic population fluctuations on diversity and total 
population sizes on evolutionary time scales are studied. In one model, species can make either 
predator-prey, mutualistic, or competitive interactions, while the other model allows only predator- 
prey interactions. When the noise in the population dynamics is strong enough, both models show 
intermittent behavior and their power spectral densities show approximate 1// fluctuations. In 
the noiseless limit, the two models have different power spectral densities. For the predator-prey 
model, l// 2 fluctuations appears, indicating random- walk like behavior, while the other model still 
shows 1/ f noise. These results indicate that stochastic population fluctuations may significantly 
affect long-time evolutionary dynamics. 
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I. INTRODUCTION 



Biological macroevolution has attracted much interest, not only among biologists, but 
also among physicists, because it is a highly nonlinear, far-from-equilibrium, complex inter- 
acting system. Recently, one of the authors proposed a series of individual-based biological 
coevolution models and studied their long-term statistical properties [2, si Q, S|- It was 
found that some properties are observed universally for various models. One such prop- 
erty is an approximate 1/f power spectral density (PSD) for the diversities (i.e., number of 
species). The models display intermittency in their time series: relatively quiet periods are 
interrupted by active periods during which the community is rearranged. During the quiet 
periods, the system is considered to be in a quasi-steady state (QSS). The 1/f fluctuations 
are so robust that one is reminded of the concept of a universality class in critical phenom- 
ena. A naturally emerging question is whether this 1/f mode is unique in evolving systems. 
In this short article, we report the existence of another mode. In that class, diversities and 
total population sizes show 1/f 2 PSDs, i.e., indicating random-walk like fluctuations. 



II. MODELS 



The models we use in this article are modified versions of individual-based coevolution 
models. The original individual-based models have a stochastic population dynamics with 
discrete, non-overlapping generations. Each species is represented by a bit-string genome 
of length L. This L-bit genome supplies a pool of 2 L possible species. At the end of 
every generation, individuals of species / give rise to F offspring before they die with the 
reproduction probability Pj. This depends on the population sizes rij{t) of all the species 
present in the community at that time. With probability (1 — Pi), the individual dies 
without offspring. The fecundity F is taken as a constant, independent of I and t for 
simplicity. In each generation, the genomes of the offspring mutate with probability \xjL 
per gene and individual. Mutation is the origin of diversity, which is necessary for the 
evolutionary process. Thus all nj(t) are updated stochastically at the same time. 

In this article, instead of individual-based population dynamics, we consider a stochastic 
difference equation to explicitly test the effect of stochasticity at the population level. In the 
case of individual-based models, the populations at t + 1 follow a binomial distribution with 
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mean ni(t)Pi and variance ni(t)Pi(l — Pj). We here approximate this process by replacing 
the binomial distribution by a gaussian distribution. Populations of continuous variables 
are updated by the following stochastic difference equation: 

ru(t + 1) = F[Pjnj(t) + Ky/m(t)Pi(l - Pi)Ut)) , (!) 

where k and are parameters giving the noise level and a random number drawn from a 
gaussian distribution with mean and standard deviation 1, respectively. The first and the 
second term in the equation represent the mean and fluctuation of the number of individuals 
in generation t + 1, respectively. When the noise coefficient k = 1, the model has the same 
mean and standard deviation as the corresponding individual-based model. The model is 
then a good approximation for an individual-based model. Populations tend to go to their 
equilibrium values, but fluctuate around them. This fluctuation may be critical, especially 
for species with tiny populations. On the other hand, when k = 0, the population dynamics 
becomes deterministic; the system converges to a fixed point and does not fluctuate after it 
reaches its fixed point. In this limit, a species can survive as long as it has a positive fixed 
point, even when its equilibrium population is quite low. 

Extinction of species is introduced by defining a threshold value, n til = 0.5. A species 
whose population becomes less than nth is considered to go extinct and is removed from 
the system. In each generation, mutations happen after the population updates with a 
probability of \xjL per gene. The number of individuals is obtained by rounding off nj(t). 

The reproduction probability is given by 

PAnAt)\ = F — ]-— t — -tt-t. (2) 

n " l + exp[-Aj({nj(t)})] 1 ; 

When Ai is large, species can almost certainly give rise to offspring; while they tend to die 
without offspring when At is small. 

n 

We use two forms of Aj, called model A and model B. Model A was introduced in 
In this model, A/ has the form 

Aj({nj(*)}) = Y,M IJ nj(t)/N tot (t) ~ N tot (t)/N Q , (3) 
j 

where Mjj, N tot , and iVo denote the interaction coefficient between I and J, the total 
population J2i n i{t), and the environmental carrying capacity, respectively. The last term 
in Eq. ([3]) limits the total population to a finite value. All the off-diagonal elements of the 
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interaction matrix Mu are randomly drawn from a uniform distribution over [—1, +1], while 
the diagonal elements Mu are set to 0. 
In model B, Aj is 

Aj(R, {nj(t)}) = -bj + V iR/Ntot(t) + £ Mijnj{t)/N^{t\ (4) 

j 

where hi and rji denote the reproduction cost and the ability to utilize the external resource 
for species /, respectively. The external resource is represented by R, which remains con- 
stant. The birth costs bi are randomly drawn from a uniform distribution over [0,1]. A 
certain proportion of species (0.05) has the ability to utilize the external resource. Thus 
7]i is positive only for these species, and it is then drawn from a uniform distribution over 
[0,1]. Other species have t]j = 0. The interaction matrix is limited to anti-symmetric 
form (Mjj = —Mji), and it is non-zero with a certain connectance probability (here, 0.1). 
Each off-diagonal Mjj is selected from a triangular distribution over [—1, 1]. The diagonal 
elements are distributed uniformly over [—1,0]. See for further details. 

III. RESULTS 

We performed Monte Carlo simulations and calcualted the PSDs of the diversity and the 
total population. We used the diversity measure known in ecology as the Shannon- Wiener 
diversity, which is defined as the exponential of the information-theoretical entropy of the 
population distribution D{t) = exp[S({rii(t)}}, where 

S({n I (t)}) = - Pi(t)ln Pl (t), (5) 

{I\ Pl (t)>0} 

with pi(t) = nj(t)/N tot (t). Since there are many kinds of species with tiny populations 
that are mostly unsuccessful mutants, this measure of diversity is useful to filter out the 
corresponding noise. 

For model A, we performed 12 independent runs for each noise level. Each simulation run 
was performed for a long period of 2 25 = 33 554 432 generations with an initializing period 
of 2 24 = 16 777216 generations. The simulation parameters were L = 13, F = 4, N = 2000, 
and fi = 0.001. Results for several noise levels are shown in Fig. [TJ It shows approximate 
1/f behavior, regardless of the noise level. However, a crossover to white noise is seen in 
the low-frequency range for strong noise, indicating the appearance of a characteristic time 
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scale. Thus, the 1/f intermittency is robust against the change in noise level in this model. 
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(a)PSD of diversities 



(b)PSD of total population sizes 



FIG. 1: Power spectral densities of diversities (a) and total population sizes (b) for model A at 
several noise levels. A line indicating 1// is also shown in both figures as a guide to the eye. 

For model B, we performed 6 independent runs for each parameter. Simulations were 
performed for 2 26 = 67108 864 generations plus a 2 24 generations initial "warm-up" period. 
The parameters were L = 18, F = 2, R = 2000, and fi = 0.0005. The results are shown in 
Fig.[2J In contrast to model A, model B exhibits PSDs that depend on the noise level. When 
the noise is strong enough, they show approximate 1/f fluctuations in a similar manner as the 
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(b)PSD of total population sizes 

FIG. 2: Power spectral densities of diversities and total population sizes with several noise levels 
for model B. Lines indicating l/f a and 1/f 2 are also shown in both figures as guides to the eye, 
with a = 1.3 and 1.5 for (a) and (b), respectively. 
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individual-based model. However, as the noise weakens, the exponents of the PSDs approach 
2; and eventually they show 1/f 2 like behavior in the noiseless limit. The exponents change 
not suddenly but gradually. A PSD with 1/f 2 behavior means that the quantities change 
like random walks, suggesting gradual changes of the species configurations. 

IV. SUMMARY AND DISCUSSION 

In this short article, we studied the effects of stochastic population fluctuations for two 
models of biological evolution. It was found that 1/f 2 PSDs of diversity and total population 
sizes may appear in model B, while 1/f PSDs are robustly found for model A. The evolving 
food webs in model B show not only 1/f, but also 1/f 2 behavior, depending on the level 
of the stochastic population fluctuations. Therefore, the species configurations for model B 
in the noiseless limit change gradually rather than intermittently, indicating the absence of 
QSS. It is desirable to analyze the noise sensitivity of model B in further detail to clarify the 
origin of the 1/f and 1/f 2 fluctuations. For example, distributions of the lifetime of species 
or durations of QSS should provide further information. 

Other forms of population dynamics including more simplistic ones and realistic 

ones jfj, have also been suggested. Comparison with these forms is a promising way to 
confirm the robustness of the results shown in this article. 
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